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A  new  empirical  atmospheric  density  model,  Jacchia-Bowman  2008,  is  developed  as  an 
improved  revision  to  the  Jacchia-Bowman  2006  model  which  is  based  on  Jacchia’s 
diffusion  equations.  Driving  solar  indices  are  computed  from  on-orbit  sensor  data  are 
used  for  the  solar  irradiances  in  the  extreme  through  far  ultraviolet,  including  x-ray  and 
Lyman-a  wavelengths.  New  exospheric  temperature  equations  are  developed  to 
represent  the  thermospheric  EUV  and  FUV  heating.  New  semiannual  density  equations 
based  on  multiple  81 -day  average  solar  indices  are  used  to  represent  the  variations  in  the 
semiannual  density  cycle  that  result  from  EUV  heating.  Geomagnetic  storm  effects  are 
modeled  using  the  Dst  index  as  the  driver  of  global  density  changes.  The  model  is 
validated  through  comparisons  with  accurate  daily  density  drag  data  previously  computed 
for  numerous  satellites  in  the  altitude  range  of  175  to  1000  km.  Model  comparisons  are 
computed  for  the  JB2008,  JB2006,  Jacchia  1970,  and  NRLMSIS  2000  models. 
Accelerometer  measurements  from  the  CHAMP  and  GRACE  satellites  are  also  used  to 
validate  the  new  geomagnetic  storm  equations. 
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I.  Introduction 


Until  development  of  the  Jacchia-Bowman  2006  (JB2006)  model1  typical  density  model  errors  on  the  order 
of  1 5%-20%  one  standard  deviation  were  recognized  for  all  empirical  models2  developed  since  the  mid 
1960s.  These  large  density  standard  deviations  correspond  to  maximum  density  errors  of  approximately 
40-60%  as  observed  in  satellite  drag  data.  There  are  two  main  reasons  for  these  consistently  large  values. 
One  is  the  result  of  not  modeling  the  semiannual  density  variation3'4  as  a  function  of  solar  activity,  and  the 
other  results  from  not  modeling  the  full  thermospheric  heating  from  solar  ultraviolet  radiation. 
Geomagnetic  storms  provide  episodic,  and  overall  smaller,  contributions  to  the  standard  deviation.  All 
models  prior  to  JB2006  have  used  the  Fi0  and  81 -day  centered  average  F,i  proxies  as  representative  of  the 

solar  ultraviolet  (UV)  heating.  However,  the  unmodeled  errors  derived  from  satellite  drag  data  all  show 
very  large  density  errors  with  approximately  27-day  periods,  representing  one  solar  rotation  cycle.  These 
errors  are  the  result  of  not  fully  modeling  the  ultraviolet  radiation  effects  on  the  thermosphere,  which  have 
a  one  solar  rotation  periodicity.  The  purpose  of  this  paper  is  to  describe  the  further  development  of  the 
Jacchia-Bowman  models  that  incorporate  new  solar  indices,  a  new  semiannual  density  model,  and  a  new 
geomagnetic  index  model. 


II.  Density  Data  Sources 

Four  different  density  data  sources  were  used  in  the  development  of  the  JB2008  model.  These  sources 
included  Air  Force  daily  density  values  from  1997  through  2007,  Air  Force  HASDM  densities  values  from 
2001  through  2005,  CHAMP  accelerometer  densities  from  2001  through  2005,  and  GRACE  accelerometer 
densities  from  2002  through  2005. 

The  daily  density  values  consist  of  very  accurate  daily  values5  obtained  from  drag  analysis  of  numerous 
satellites  with  perigee  altitudes  of  175  km  to  1000  km.  Daily  temperature  corrections  to  the  Air  Force 
modified  Jacchia  1970  atmospheric  model6,7  were  obtained  from  the  satellite  data  throughout  the  period 
1997  through  2007.  Approximately  210,000  daily  temperature  values  were  computed  using  a  special 
energy  dissipation  rate  (EDR)  method5,  where  radar  and  optical  observations  are  fit  with  special  orbit 
perturbations.  For  each  satellite  tracked  approximately  100,000  radar  and  optical  observations  were 
available  for  the  special  perturbation  orbit  fitting.  A  differential  orbit  correction  program  was  used  to  fit 
the  observations  to  obtain  the  standard  6  Keplerian  elements  plus  the  ballistic  coefficient.  “True”  ballistic 
coefficients8  were  then  used  with  the  observed  daily  temperature  corrections  to  obtain  daily  density  values. 
The  daily  density  computation  was  validated  by  comparing  historical  daily  density  values  computed  for  the 
last  30  years  for  over  30  satellites.  The  accuracy  of  the  density  values  was  determined  from  comparisons  of 
geographically  overlapping  perigee  location  data,  with  over  8500  pairs  of  density  values  used  in  the 
comparisons.  The  density  errors  were  found  to  be  less  than  4%  overall,  with  errors  on  the  order  of  2%  for 
values  covering  the  latest  solar  maximum. 

Density  values  were  also  obtained  from  the  HASDM9,10  model.  The  Air  Force  Space  Command’s  High 
Accuracy  Satellite  Drag  Model  (HASDM)  processes  drag  information  from  the  trajectories  of  75  to  80 
inactive  payloads  and  debris  (calibration  satellites)  to  solve  for  a  dynamically  changing  global  correction  to 
the  thermosphere.  This  correction  covers  the  altitude  range  of  -200  to  800  km.  Satellite  tracking 
observations  (azimuth,  elevation,  range,  and  range  rate)  from  the  Air  Force  Space  Surveillance  Network 
(SSN)  are  processed  directly  to  derive  the  neutral  atmospheric  density.  Thermospheric  density  correction 
parameters  are  computed  along  with  the  trajectory  states  of  the  calibration  satellites  in  a  single  estimation 
process,  known  as  the  Dynamic  Calibration  Atmosphere  (DCA).  This  is  a  weighted  least  squares  differential 
correction  across  all  calibration  satellites  that  simultaneously  solve  for  global  density  corrections  and  a  state 
vector  for  each  calibration  satellite.  DCA  uses  the  Jacchia  1970  thermosphere  as  its  base  model.  DCA 
estimates  13  global  temperature  correction  parameters  every  3  hours,  while  the  state  vector  of  each 
calibration  satellite  is  estimated  for  a  2-day  fit  span.  The  13  global  temperature  correction  parameters 
produce  a  spatial  resolution  of  approximately  80  degree  geographical  blocks.  The  HASDM  and  DCA 
programs  have  been  previously  well  documented9'10.  For  the  JB2008  model  development  densities  were 
computed  every  10  seconds  along  the  CHAMP  and  GRACE  orbits  using  the  HASDM  temperature 
coefficients  obtained  for  the  2001  through  2005  time  period. 
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Another  density  data  source  came  from  the  CHAMP  (CHAllenging  Minisatellite  Payload)  satellite,  a 
German  small  satellite  mission  for  geoscientific  and  atmospheric  research  and  applications,  managed  by 
GFZ11,  Potsdam.  CHAMP  was  launched  on  July  15,  2000  into  an  almost  circular,  near-polar  orbit 
(inclination  87.2°)  with  an  initial  average  altitude  of  450  km.  CHAMP  carries  a  very  sensitive  STAR 
accelerometer,  the  data  of  which  can  be  used  to  derive  neutral  densities12.  Densities  every  10  seconds  were 
available  for  the  2001  through  2005  time  period. 


A  fourth  density  data  source  used  in  this  model  development  came  from  the  GRACE  satellite  mission 
(Gravity  Recovery  and  Climate  Experiment)13,  the  mission  objective  being  to  map  the  global  gravity  field 
with  unprecedented  accuracy.  GRACE  is  a  twin  satellite  configuration,  which  was  launched  on  March  1 7, 
2002  into  an  almost  circular,  near-polar  orbit  (inclination  89.0°)  with  an  initial  altitude  of  500  km.  GRACE 
carries  extremely  sensitive  SuperSTAR  accelerometers  which  are  an  order  of  magnitude  more  precise  than 
STAR.  Densities  every  5  seconds  were  available14  for  the  2002  through  2005  time  period. 

III.  Global  Nighttime  Minimum  Exospheric  Temperature 
A.  Temperature  Description 

The  variations  in  the  ultraviolet  solar  radiation  that  heats  the  earth's  thermosphere  consists  of  two 
components,  one  related  to  solar  rotational  modulation  of  active  region  emission,  and  the  other  long-term 
evolution  of  the  main  solar  magnetic  field6  7.  The  passage  of  active  regions  across  the  disk  during  a  solar 
rotation  period  produces  irradiance  variations  of  approximately  27  days,  while  the  main  solar  magnetic 
field  evolution  produces  irradiance  variations  over  approximately  11  years.  The  10.7-cm  solar  flux,  F10, 
has  in  the  past  been  used  to  represent  these  effects.  However,  new  solar  indices  have  been  recently15  used 
to  compute  better  density  variation  correlations  with  ultraviolet  radiation  covering  the  entire  Far  UV  as  well 
as  the  EUV  wavelengths. 

In  determining  a  new  global  nighttime  minimum  exospheric  temperature  Tc  equation  with  the  new  solar 
indices,  the  density  values  were  converted  into  daily  Tc  temperature  values  using  the  Jacchia  70  empirical 
atmospheric  density  model.  To  obtain  accurate  Tc  values  the  large  semiannual  density  variations  had  to  be 
correctly  modeled.  A  major  density  variation,  aside  from  the  1 1-year  and  27-day  solar  heating  effect,  is  the 
semiannual  change.  This  can  be  as  large  as  250%  from  a  July  minimum  to  an  October  maximum  during 
solar  maximum  years,  and  as  small  as  60%  from  July  to  October  during  solar  minimum  years  (at  600  km)3. 
The  semiannual  variation  was  computed  on  a  yearly  basis  from  the  previously  derived  density  data. 
Jacchia’s  70  semiannual  density  model  equation  was  then  replaced  using  these  observed  semiannual  yearly 
variations.  A  smaller  correction  to  Jacchia’s  model  was  also  made  for  the  observed  errors  in  the  latitude 
and  local  solar  time  density  variations.  From  these  different  model  corrections  an  accurate  Tc  value,  due 
almost  entirely  to  solar  heating,  was  obtained. 

The  solar  UV  absorption  spectrum  in  the  thermosphere  was  analyzed  to  determine  the  new  solar  indices 
required  for  the  new  temperature  equation  development.  The  solar  index  Fi0  is  really  a  proxy  index 
because  it  is  measured  at  a  10.7-cm  wavelength,  which  is  not  a  direct  measure  of  any  ultraviolet  radiation 
and  is  not  absorbed  by  the  atmosphere.  Direct  ultraviolet  heating  indices  were  recently  developed  that 
represent  the  extreme  (EUV),  far  (FUV),  and  mid  (MUV)  solar  UV  radiation.  Previous  analyses15 
suggested  that  EUV  and  FUV  indices26  were  required  to  capture  most  of  the  thermospheric  heating,  and  an 
additional  improvement  could  be  obtained  by  using  an  index  representing  UV  energy  absorption  at  lower 
thermospheric  altitudes  than  by  using  previous  EUV  and  FUV  indices.  The  daily  indices  selected  for  this 
model  development  include  F[0,  Si0,  M10,  and  Yi0. 

F i0:  The  10.7-cm  solar  radio  flux,  F10,  is  produced  daily  by  the  Canadian  National  Research  Council’s 
Herzberg  Institute  of  Astrophysics  at  its  ground-based  Dominion  Radio  Astrophysical  Observatory  located 
in  Penticton,  British  Columbia.  The  physical  units  of  Fj0  are  W  m"2  Hz"1  and  more  conveniently  expressed 
in  solar  flux  units  (1  sfu  =  lxlO"22  W  m"2  Hz"1).  A  running  81-day  centered  smoothed  set  of  values  using 
the  moving  boxcar  method  was  created,  and  these  data  are  referred  to  as  F10 .  Linear  regression  with  daily 
F  io  has  been  used  to  scale  and  report  all  other  solar  indices  in  units  of  sfu. 
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Sio:  The  NASA/ESA  Solar  and  Heliospheric  Observatory  (SOHO)  research  satellite  uses  the  Solar 
Extreme-ultraviolet  Monitor  (SEM)  instrument  that  has  been  measuring  the  26-34  nm  solar  EUV  emission 
since  launch  in  December  1995.  This  integrated  26-34  nm  emission,  which  is  also  measured  by  the  post- 
GOES  12  operational  satellites,  has  been  normalized  and  converted  to  sfu  through  linear  regression  with 
Fio,  producing  the  new  index  Si0.  The  broadband  (wavelength  integrated)  SEM  26-34  nm  irradiances  are 
EUV  line  emissions  dominated  by  the  chromospheric  He  II  line  at  30.4  nm  with  contributions  from  other 
chromospheric  and  coronal  lines.  This  energy  principally  comes  from  solar  active  regions.  A  running  81- 
day  centered  smoothed  set  of  values  using  the  moving  boxcar  method  was  created,  and  these  data  are 
referred  to  as  S10 . 

Mi0:  The  NOAA  series  of  operational  satellites  host  the  Solar  Backscatter  Ultraviolet  (SBUV) 
spectrometer  that  has  the  objective  of  monitoring  ozone  in  the  Earth’s  lower  atmosphere.  In  its  discrete 
operating  mode,  a  diffuser  screen  is  placed  in  front  of  the  instillment’s  aperture  in  order  to  scatter  solar 
MUV  radiation  near  280  nm  into  the  instrument.  This  solar  spectral  region  contains  both  photospheric 
continuum  and  chromospheric  line  emissions.  The  chromospheric  Mg  II  h  and  k  lines  at  279.56  and  280.27 
nm,  respectively,  and  the  weakly  varying  photospheric  wings  (or  continuum  longward  and  shortward  of  the 
core  line  emission),  are  operationally  observed  by  the  instrument.  The  Mg  II  core-to-wing  ratio  (cwr)  is 
calculated  between  the  variable  lines  and  nearly  non-varying  wings.  The  result  is  a  measure  of 
chromospheric  and  some  photospheric  solar  active  region  activity  independent  of  instrument  sensitivity 
change  through  time,  and  is  referred  to  as  the  Mg  II  cwr.  The  Mg  II  cwr  have  been  used  in  a  linear 
regression  with  F10  to  derive  the  M10  index  in  sfu  units.  A  running  81-day  centered  smoothed  set  of  values 
using  the  moving  boxcar  method  was  created,  and  these  data  are  referred  to  as  M  . 

Y10:  The  operational  GOES  X-ray  Spectrometer  (XRS)  instrument  provides  the  0. 1-0.8  nm  solar  X-ray 
emission.  X-rays  in  the  0. 1-0.8  nm  range  come  from  the  cool  and  hot  corona  and  are  typically  a 
combination  of  both  very  bright  solar  active  region  background  that  varies  slowly  (days  to  months)  plus 
flares  that  vary  rapidly  (minutes  to  hours),  respectively.  The  photons  arriving  at  Earth  are  primarily 
absorbed  in  the  mesosphere  and  lower  thermosphere  (80-90  km)  by  molecular  oxygen  and  nitrogen  where 
they  ionize  those  neutral  constituents  to  create  the  ionospheric  D-region.  An  index  of  the  solar  X-ray  active 
region  background,  without  the  flare  component,  has  been  developed.  This  is  called  the  X]0  index27.  The 
0. 1-0.8  nm  X-rays  are  a  major  energy  source  in  these  atmospheric  regions  during  high  solar  activity  but 
relinquish  their  dominance  to  the  competing  hydrogen  (H)  Lyman-a  emission  during  moderate  and  low 
solar  activity.  Lyman-a  is  also  deposited  in  the  same  atmospheric  regions,  is  created  in  the  solar  upper 
chromosphere  and  transition  region,  and  demarcates  the  EUV  from  the  FUV  spectral  regions28.  It  is  formed 
primarily  in  solar  active  regions,  plage,  and  network;  the  photons,  arriving  at  Earth,  are  absorbed  in  the 
mesosphere  and  lower  thermosphere  where  they  dissociate  nitric  oxide  (NO)  and  participate  in  water  (H20) 
chemistry.  Lyman-a  has  been  observed  by  the  SOLSTICE  instrument  on  the  UARS  and  SORCE  NASA 
research  satellites  as  well  as  by  the  SEE  instrument  on  NASA  TIMED  research  satellite.  Since  these  two 
solar  emissions  are  competing  drivers  to  the  mesosphere  and  lower  thermosphere,  we  have  developed  a 
mixed  solar  index  Yi0.  It  is  weighted  to  represent  mostly  X10  during  solar  maximum  and  to  represent 
mostly  Lyman-a  during  moderate  and  low  solar  activity.  The  independent,  normalized  ij0  is  used  as  the 
weighting  function  and  multiplied  with  the  Xi0  and  Lyman-a  as  fractions  to  their  solar  maximum  values. 

B.  Tc  Temperature  Equation 

Previous  analyses16'17  of  different  density  model  errors  have  shown  that  using  the  Fj0  index  to  capture  the 

11-year  solar  cycle  variation  does  not  fully  represent  the  entire  thermospheric  heating,  especially  during 
solar  minimum  conditions.  It  has  been  shown  that  real  density-to-model  ratios  have  drops  of  30-40%  at 
solar  minimum.  The  F(.  index  has  long  been  known  to  “flatten-ouf  ’  around  solar  minimum,  while  the  real 

EUV  heating  continues  to  show  variability.  However,  a  previous  analysis15  demonstrated  that  the  Fj.  index 

was  still  better  at  representing  the  full  1 1-year  cycle  changes  than  either  the  S  or  M  index.  Therefore,  it 

was  decided  to  use  the  Fjn  index  for  the  great  majority  of  the  time,  but  supplementing  it  with  the  EUV  S10 

during  solar  minimum  times.  With  this  approach  a  new  11-year  solar  index  was  developed  with  the 
following  weighting  scheme: 
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(1) 


Fs  =  4  x  WT  +  S10  x  (1-WT)  where  WT  =  ^J240) 

With  this  new  index  the  solution  of  the  best  nighttime  minimum  exospheric  Tc  equation  was  obtained  using 
numerous  satellites  for  the  years  from  1997  through  2007  when  all  new  solar  indices  were  available.  The 
resulting  equation  is: 

Tc  =  392.4  +  3.227  Fs  +  0.298  AF10  +  2.259  AS10  +  0.3 12  AM10  +  0. 1 78  AY10  (2) 

The  delta  values  (AFjo,  ASjo,  AMjo,  AY10)  represent  the  difference  of  the  daily  and  81-day  centered  average 
value  of  each  index.  The  81-day  (3-solar  rotation  period)  centered  value  was  previously15  determined  to  be 
the  best  long  term  average  to  use.  In  the  solution  the  2007  solar  minimum  data  was  heavily  weighted  to 
help  better  define  the  density  variations  during  solar  minimum  times.  To  avoid  increases  in  Tc  due  to 
geomagnetic  storms  all  daily  data  with  the  geomagnetic  index  %>=  30  were  rejected.  This  meant  that  if  a 
solar  index  required  a  lag  time  of  5  days,  each  of  the  5  days  prior  to  the  current  time  had  to  have  ap  <  30  for 
the  current  daily  density  data  to  be  used. 

It  was  determined  that  a  lag  time  of  1  day  was  the  best  to  use  for  the  Fi0  and  Sio  indices.  For  using  the  Mi0 
index  an  analysis  determined  that  the  best  (least  squares  minimum)  lag  time  was  2  days,  and  for  Y10  a  best 
lag  time  of  5  days  was  obtained. 

Initially  for  the  JB2006  model,  which  did  not  use  Y]0,  the  lag  time  for  M10  was  determined  to  be  5  days. 
The  M10  index  was  previously  accounting  for  the  longer  lag  times  in  the  lower  thermosphere.  Flowever, 
with  the  addition  of  the  low  altitude  Y 10  index  the  M10  lag  time  became  shorter,  and  the  low  altitude  longer 
absorption  lag  time  was  captured  by  Y 10  combining  absorption  of  X-Rays  and  Lyman-a  at  altitudes  around 
80-90  km. 

In-order  to  evaluate  the  new  Tc  equation  the  “observed”  density-to-model  ratios  were  computed  for  both 
the  JB2006  and  new  JB2008  models,  the  Jacchia  70  model,  and  the  NRLMSIS18  model.  The  new  JB2008 
semiannual  equations,  discussed  in  the  following  sections,  were  used  in  the  JB2008  evaluation.  The 
“observed”  densities  were  obtained  by  using  the  computed  3 -hour  spherical  harmonic  FIASDM  temperature 
correction  coefficients,  and  computing  density  values  at  10  minute  steps  along  the  CF1AMP  reference  orbits 
obtained  for  2001  through  2007.  These  FlASDM-to-Model  ratios  were  then  binned  by  Fj.  and  plotted  in 

Figure  1.  It  can  be  readily  seen  that  all  the  previous  models  using  just  1^0  for  the  11 -year  cycle  variations 

show  a  significant  decrease  in  the  ratios  at  solar  minimum  conditions.  The  JB2008  model  does  much 
better  at  representing  the  solar  minimum  density  decrease,  although  it  still  does  not  completely  capture  the 
density  variation.  Figure  2  shows  the  density  model  standard  deviations  binned  again  by  Ij(1 .  The  much 

larger  sigma  at  solar  minimum  (very  low  Ij() )  are  a  direct  result  of  the  model  ratio  errors  at  low  Fy  .  The 

new  JB2008  Tc  equation  is  a  significant  improvement  over  all  other  models  in  representing  the  solar 
thermospheric  heating. 
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HASDM  /  MODEL  Density  Ratio  vs  F10B  During  2001-2007  (  400  km  ) 


Figure  1.  HASDM-to-Model  density  ratios  at  400km  altitude  as  a  function  of  Fi:j  (F10B). 


Density  STD  vs  FI  0B  During  2001-2007  (  400  km  ) 


Figure  2.  Density  percentage  errors  (1  standard  deviation)  from  model  density  values  at  400  km  altitude 
compared  to  HASDM  density  values. 
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IV.  Semiannual  Density  Variation 


The  semiannual  density  variation  was  first  discovered  in  1961 19.  Paetzold  and  Zschorner  observed  a  global 
density  variation  from  analysis  of  satellite  drag  data,  which  showed  a  6-month  periodicity  maximum 
occurring  in  April  and  October,  and  minimum  occurring  in  January  and  July.  Many  authors,  such  as  King- 
Hele20  and  Jacchia6’7,  analyzed  the  semiannual  effect  from  satellite  drag  during  the  1960s  and  early  1970s. 
They  found  that  the  semiannual  variation  was  a  worldwide  effect  with  the  times  of  the  yearly  maximum 
and  minimum  occurring  independent  of  height.  However,  the  semiannual  period  was  found  to  be  only 
approximate,  as  the  times  of  occurrence  of  the  minima  and  maxima  seemed  to  vary  from  year  to  year. 
Generally  the  October  maximum  exceeded  that  in  April  and  the  July  minimum  was  deeper  than  that  in 
January.  The  main  driving  mechanism  for  the  observed  variability  in  the  semiannual  effect  remained  a 
mystery.  Jacchia  first  modeled  the  effect  as  a  temperature  variation  which  included  a  function  of  the  81- 

day  solar  flux  Fj0  index.  However,  he  soon  discovered  difficulties  with  the  temperature  model,  and 
eventually  modeled  the  semiannual  variation  as  a  density  variation.  He  also  dropped  the  F|(j  dependence, 

suggesting  that  he  did  not  have  enough  data  to  support  this  solar  flux  relationship.  He  found  that  the 
amplitude  of  the  semiannual  density  variation  was  strongly  height-dependent  and  variable  from  year  to 
year.  However,  he  could  not  show  a  definitive  correlation  of  the  variation  with  solar  activity. 

A.  Semiannual  Density  Variation  Function 

Jacchia  obtained  the  following  equations  from  analysis  of  12  years  of  satellite  drag  data.  He  represented 
the  semiannual  density  variation  in  the  form: 

asa  iogio  P  =  F(z)  G(t)  (3) 

G(t)  represents  the  average  density  variation  as  a  function  of  time  in  which  the  amplitude  (i.e.  the 
difference  in  log10  density  between  the  principal  minimum  in  July  and  the  principle  maximum  in  October) 
is  normalized  to  1,  and  F(z)  is  the  relation  between  the  amplitude  and  the  height  z. 

From  previous  analysis3  it  was  determined  that  a  Fourier  series  could  accurately  represent  Jacchia’s  G(t) 
equation  structure.  A  9-term  coefficient  series,  including  frequencies  up  to  4  cycles  per  year,  was  sufficient 
to  capture  all  the  yearly  variability  in  G(t)  that  had  been  previously  observed  by  Jacchia.  It  was  also 
determined  that  a  simplified  quadratic  polynomial  equation  in  z  could  sufficiently  capture  Jacchia’s  F(z) 
equation  and  not  lose  any  fidelity  in  the  observed  F(z)  values. 

B.  Semiannual  F(z)  Height  Function 

For  the  Jacchia-Bowman  model  developments  the  amplitude,  F(z),  of  the  semiannual  variation  was 
determined3  4  on  a  year-by-year  and  satellite-by-satellite  basis.  The  smoothed  density  difference  data  was 
fit  each  year  for  each  satellite  using  a  9  term  Fourier  series.  The  F(z)  value  was  then  computed  from  each 
fit  as  the  difference  between  the  minimum  and  maximum  values  for  the  year. 

Figure  3  shows  the  results  of  three  different  years  of  data,  along  with  the  plot  of  Jacchia’s  F(z)  equation. 
For  each  year,  the  F(z)  values  were  fit  with  a  quadratic  polynomial  in  height.  The  smoothed  curves  shown 
in  Figure  3  represent  the  least  squares  quadratic  fit  obtained  for  three  different  years.  The  F(z) 
Alog10p  data  for  all  satellites  are  very  consistent  within  each  year.  The  most  notable  feature  in  Figure  3  is 

the  very  large  difference  in  maximum  amplitude  among  the  years  displayed.  The  2002  data  shows  a 
maximum  density  variation  of  250%  near  800km,  while  the  1993  data  shows  only  a  60%  maximum 
variation.  Jacchia’s  F(z)  function  only  gives  a  constant  130%  maximum  variation  for  all  years. 

Previous  development  of  the  JB2006  model  showed  that  solar  EUV  and  FUV  heating  played  an  important 
part  in  thermospheric  density  variations.  Bowman4  extended  the  previous  semiannual  work3  to  include 
additional  solar  EUV  indices  in  an  attempt  to  capture  the  remaining  semiannual  variations  not  modeled  by 
the  JB2006  model. 
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Figure  3.  The  amplitude  function  F(z)  for  three  different  years  (1990,  1993,  2002),  with  semiannual 
amplitudes  plotted  for  each  satellite  for  each  year.  The  constant  F(z)  function  from  Jacchia  is  also  plotted. 


Roble21  computed  the  thermospheric  temperature  response  to  solar  EUV  heating  using  his  coupled 
thermosphere  and  ionosphere  global  average  model.  Fie  found  that  removing  the  Fie  II  30.4  nm  emission 
produced  the  largest  (by  a  factor  of  2)  temperature  change.  Therefore,  it  was  very  important  to  select  an 
EUV  index  that  captured  the  emission  of  this  He  II  irradiance  line.  These  results  together  with  previous 
analysis15  of  thermospheric  response  to  new  solar  indices  suggested  a  new  set  of  solar  indices  to  use  for  the 

semiannual  variation.  New  81 -day  centered  S  and  M  indices  were  computed  for  use  along  with  the 
previous  F.  index.  Previous  work4  determined  the  new  solar  index  for  F(z)  to  be 

Fsmj  -  1-00  Fj  -  0.70  S3  -  0.04  MJ  (4) 

where  the  Fj ,  Sj ,  and  M;  indices  represent  the  July  averages  of  the  Fj0 ,  S10 ,  and  M10  indices 
respectively.  This  FSMJ  index  was  then  used  to  determine  which  terms  were  significant  in  defining  a  new 

F(z)  equation.  The  resultant  new  F(z)  equation,  with  z  =  height/1000,  using  the  new  index  was  determined 
to  be 


F(z)  -  B,  +  B2FSMJ  +  B3zFSMJ  +  B4z2FSMJ  +  B5zFgMJ  (5) 

Table  1  lists  the  resulting  B  coefficient  values  with  their  standard  deviations  obtained  from  using  Equation 
(4)  for  the  solar  index  required  in  Equation  (5).  The  standard  deviations  of  all  the  coefficients  are  an  order 
of  magnitude  less  than  the  coefficient  values,  indicating  that  all  five  coefficients  have  been  well 
determined. 

Equation  (5)  using  FSMJ  represents  a  global  equation  in  F(z)  using  data  from  yearly  semiannual  amplitudes 
observed  from  1997  through  2006.  For  incorporation  into  JB2008  the  81-day  centered  July  average  values 
are  replaced  by  daily  81-day  centered  values  of  IJ0 ,  SI0 ,  and  M]0 .  This  is  an  approximation  to  the  best  fit 

equation.  Using  the  daily  81 -day  centered  values  in  Equations  (4)  and  (5)  result  in  an  increase  in  the 
density  error  standard  deviation  of  less  than  1%. 
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Coef 

Term 

Value 

STD 

Bi 

1 

2.69E-01 

1.84E-02 

b2 

F 

SM  J 

-1.18E-02 

6.56E-04 

b3 

z  F 

A  ±SMJ 

2.78E-02 

1.92E-03 

b4 

2  p 

7  r 
^  XSMJ 

-2.78E-02 

1.20E-03 

b5 

c2 

zFsmj 

3.47E-04 

3.51  E-05 

Table  1.  F(z)  coefficient  values  with  standard  deviations  (STD)  from  best  fit  results. 


C.  Semiannual  G(t)  Yearly  Periodic  Function 

The  yearly  observed  G(t)  function,  as  previously  discussed,  consists  of  a  Fourier  series  with  9  coefficients 
representing  a  quadannual  variation.  28-day  smoothed  density  difference  data  for  each  satellite  was  fitted 
with  this  Fourier  series  for  each  year.  The  density  difference  data  is  the  accurate  observed  daily  density 
values  minus  the  Jacchia  values  without  Jacchia’s  semiannual  variation.  The  G(t)  function  was  then 
obtained  by  normalizing  to  a  value  of  1.0  the  difference  between  the  minimum  and  maximum  values  for  the 
year.  The  F(z)  value  for  each  satellite  by  year  was  used  for  the  normalization.  Figure  4  shows  the  results 
obtained  for  the  year  1990  for  the  majority  of  the  satellites.  Note  the  tight  consistency  of  the  curves  for  all 
heights,  covering  over  800  km  in  altitude,  which  demonstrates  the  validity  of  using  one  G(t)  function  per 
year  to  represent  the  yearly  semiannual  phase  for  all  altitudes.  This  tight  consistency  of  the  G(t)  phase  for 
all  satellites  also  indicates  that  there  is  no  significant  latitude  or  local  solar  time  effects  with  the  semiannual 
density  variation.  This  conclusion  can  be  made  because  the  majority  of  the  satellites  have  moderate  to  high 
eccentricity  orbits.  This  means  that  the  great  majority  of  the  density  sampling  on  each  revolution  occurs 
very  close  to  the  perigee  location,  and  the  daily  density  values  computed  from  the  orbit  decays  can  be 
assigned  to  the  argument  of  perigee  latitude  and  local  solar  time,  which  is  different  for  each  satellite.  The 
precession  of  the  argument  of  perigee  can  be  very  slow  (from  zero  to  a  few  degrees  per  day),  so  if  there  is  a 
latitude  or  local  solar  time  semiannual  effect  the  G(t)  phase  curves  in  Figure  4  should  show  significant 
differences  because  of  the  random  nature  of  the  argument  of  perigee  locations.  This  is  definitely  not 
observed  when  comparing  all  of  the  individual  satellite  G(t)  phase  curves. 

The  next  step  in  the  study  was  to  fit  a  yearly  9-term  G(t)  function  for  each  year  using  the  data  for  all  the 
satellites  for  the  year.  Figure  4  also  shows  the  yearly  fit  G(t)  value  for  the  year  1990.  A  small  standard 
deviation  was  obtained  for  every  year’s  fit,  especially  during  solar  maximum  years.  Figure  5  shows  the 
yearly  G(t)  fits  for  1999  through  2001,  again  showing  the  consistency  of  the  semiannual  phase  at  all 
altitudes  for  a  given  year.  Also,  it  is  readily  apparent  that  the  series  changes  dramatically  from  year  to  year. 
It  was  determined  that  during  solar  maximum  the  July  minimum  date  can  vary  by  as  much  as  80  days. 
During  solar  minimum  the  semiannual  July  minimum  time  variation  is  much  smaller  and  appears  to  be 
flattened  out  in  time. 

As  was  done  for  the  F(z)  analysis  it  was  decided  to  combine  the  new  81 -day  average  indices  in  a  linear 
function  since  each  index  is  expressed  in  terms  of  Fj0  units  and  this  approach  worked  very  well  for  the  F(z) 
analysis.  A  new  solar  index,  representing  long  term  EUV  and  FUV  heating,  was  determined  to  be 

Fsm=  1.00  Fw  -  0.75  S10  -  0.37  M10  (6) 

It  was  decided  to  start  out  using  only  annual  and  semiannual  terms,  instead  of  the  JB2006  quadannual  terms 
previously  used,  to  try  to  represent  the  yearly  semiannual  phase  variations.  The  yearly  observed  values  had 
been  fit  with  terms  up  to  quadannual,  but  it  was  hoped  that  only  terms  up  to  semiannual  needed  to  be 
included  for  a  global  model.  The  resulting  equation  was 

G(t)  =  Cj  +  C2  sin(&>)  +  C3  cos(<a)  +  C4  sin(2&>)  +  C5  cos(2&>)  ^ 

+  Fsm  {  C6  +  C7  sin(ry)  +  C8  cos(&>)  +  C9  sin(2&>)  +  C10  cos(2&>)} 
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Figure  4.  The  individual  satellite  G(t)  fits  are  plotted  for  1990.  The  Jacchia  model  and  yearly  fit 
model  are  also  shown. 


The  coefficients  in  Equation  (6)  are  better  defined  than  those  for  the  F(z)  index  function  specified  by 
Equation  (4).  This  is  because  density  (G(t))  data  and  FSM  values  were  available  throughout  the  entire  year 
as  opposed  to  using  one  July  averaged  value  per  year  to  derive  Equation  (4). 

Table  2  lists  the  resulting  C  coefficient  values  with  their  standard  deviations  obtained  from  using  Equation 
(6)  for  the  solar  index  used  in  Equation  (7).  The  standard  deviations  of  the  coefficients  are  all  an  order  of 
magnitude  smaller  than  the  coefficient  values  except  for  the  C7  and  Cg  FSM  annual  terms,  indicating  a  well 
determined  set  of  coefficients. 


Coef 

Term 

Value 

STD 

Cl 

1 

-3.63E-01 

6.33E-03 

c2 

sin(ry ) 

8.51  E-02 

9.23E-03 

C3 

cos  (co) 

2.40E-01 

8.60E-03 

C4 

sin  (2  co) 

-1.90E-01 

8.61  E-03 

Cs 

cos(2  co) 

-2.55E-01 

8.79E-03 

Ce 

f" 

A  SM 

-1.79E-02 

3.63E-04 

C7 

FSMsin(®) 

5.65E-04 

5.39E-04 

Ca 

FSMc°s(o) 

-6.41E-04 

4.77E-04 

C9 

Fsm  sin(2ffi) 

-3.42E-03 

4.91  E-04 

ClO 

4,008(2®) 

-1.25E-03 

5.07E-04 

Table  2.  G(t)  coefficient  values  with  standard  deviations  (STD)  from  the  best  fit  results. 
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Figure  5.  The  individual  satellite  fits  for  3  different  years  is  shown.  The  Year  G(t)  Model  is 
highlighted.  Each  set  of  curves  for  1999  and  2001  has  been  offset  by  +1.00  and  -1.00 
respectively  in  G(t)  for  clarity.  The  JB2006  and  new  JB2008  model  curves  are  also  displayed. 


The  results  of  the  new  global  model  from  Equations  (6)  and  (7)  are  plotted  in  Figure  5  as  the  FSMB  model. 
Also  plotted  are  the  yearly  observed  values  for  each  year,  and  the  original  JB2006  global  model  values. 

The  10-term  new  model  results  are  impressive.  Even  with  only  annual  and  semiannual  terms  the  new 
model  accounts  almost  completely  for  the  July  minimum  phase  shifting  which  could  not  be  captured  in  the 
F()  global  model  using  even  quadannual  terms.  This  clearly  demonstrates  that  the  large  majority  of  the 
variations  observed  in  the  semiannual  density  variation  can  be  attributed  to  direct  solar  heating  responses. 


V.  Geomagnetic  Storm  Modeling 


A.  Dst  Index  Description 

The  Disturbance  Storm  Time  (Dst)  index  is  primarily  used  to  indicate  the  strength  of  the  storm-time  ring 
current  in  the  inner  magnetosphere.  During  the  main  phase  of  magnetic  storms,  the  ring  current  becomes 
highly  energized  and  produces  a  southward-directed  magnetic  field  perturbation  at  low  latitudes  on  the 
Earth’s  surface.  This  is  opposite  to  the  normal  northward-directed  main  field.  The  index  is  determined 
from  hourly  measurements  of  the  magnetic  field  made  at  four  points  around  the  Earth’s  equator  and  is 
available  from  World  Data  Center  (A)  at  Kyoto2". 

Most  magnetic  storms  begin  with  sharp  rises  in  Dst,  called  the  storm  sudden  commencement,  in  response  to 
increased  solar  wind  pressure.  Following  a  southward  turning  of  the  interplanetary  magnetic  field,  Dst 
decreases  as  ring  current  energy  increases  during  the  storm’s  main  phase.  During  the  recovery  phase  the 
ring  current  energy  decreases  and  Dst  increases  until  the  storm’s  end.  Traces  of  Dst  show  a  transition  from 
the  early  to  late  recovery  phase  characterized  by  significant  changes  in  slope  as  the  distribution  of  the  ring 
current  becomes  symmetric  in  local-time.  However,  a  significant  fraction  of  magnetic  storms  manifest 
more  complex  structuring,  with  multiple  main  and  partial  recovery  phases.  Figure  6  for  an  example  of  the 
Dst  events  during  a  complex  storm. 


11 

American  Institute  of  Aeronautics  and  Astronautics 


Use  of  Dst  as  a  parameter  of  the  energy  deposited  in  the  thermosphere  during  magnetic  storms  is  more 
accurate  than  the  use  of  the  ap  index.  The  3 -hour  ap  is  an  indicator  of  general  magnetic  activity  over  the 
Earth  and  responds  primarily  to  currents  flowing  in  the  ionosphere  and  only  secondarily  to  magnetospheric 
variations.  The  ap  index  is  determined  by  observatories  at  high  latitudes  which  can  be  blind  to  energy 
input  during  large  storms23  and  thus  underestimate  the  effects  of  storms  on  the  thermosphere. 

As  described  below  the  thermosphere  acts  during  storm  periods  as  a  driven-but-dissipative  system  whose 
dynamics  is  represented  by  a  differential  equation,  with  the  changes  in  exospheric  temperature  change 
given  as  a  function  of  Dst.  To  determine  the  exospheric  temperature,  and  thereby  the  thermospheric 
density  distribution  at  any  time  in  a  storm,  it  is  necessary  to  integrate  the  differential  equation  for  dTc 
starting  at  the  storm  commencement  and  proceeding  throughout  the  entire  storm  period.  Therefore,  it  is 
necessary  to  recognize  where  Dst  measurements  come  in  a  particular  storm’s  development. 

An  algorithm  for  determine  the  storm  events  was  developed  for  locating  in  time  the  start,  Dst  minimum, 
recovery  slope  change,  and  final  end  of  the  storm.  For  practical  reasons  we  define  a  magnetic  disturbance 
as  a  storm  only  if  the  minimum  Dst  <  -75  nT.  We  selected  this  value  because  disturbances  with  minimum 
Dst  >  -75  nT  often  lack  identifiable  storm  profiles.  Once  the  starting  point  of  the  storm  is  determined  the 
algorithm  steps  forward  in  time  until  the  minimum  Dst  value  is  obtained.  This  is  defined  as  the  end  of  the 
storm  main  phase.  Because  individual  Dst  traces  may  exhibit  several  local  minima  before  reaching  the 
deepest  minimum,  the  algorithm  specifies  the  real  storm  minimum  point.  Once  the  minimum  is  identified 
the  algorithm  continues  stepping  forward  through  the  recovery  phase  until  a  major  slope  change  is  detected. 
From  this  point  to  the  end  of  the  storm  the  Dst  slope  is  relatively  shallow.  It  has  been  found  that  Dst  takes 
much  longer  to  recover  than  does  the  thermosphere24.  To  determine  a  “real”  density  recovery  time  more 
than  80  storms  were  analyzed.  A  linear  fit  of  storm  duration  verse  storm  magnitude  was  obtained  to  give 
an  equation  for  the  approximate  end  time  of  the  storm.  The  algorithm  determines  if  the  storm  ends  before 
this  by  examining  when  the  Dst  values  are  above  the  -75  nT  limit.  The  lesser  in  time  of  the  Dst  limit  or 
linear  fit  time  is  used  for  the  end  time.  For  complex  storms  (a  second  disturbance  starts  before  the  previous 
one  ends)  the  algorithm  determines  the  start,  minimum,  recovery  slope  change,  and  end  point  events  of 
each  storm.  For  a  multiple  storm  the  starting  time  of  the  second  storm  will  be  at  the  same  time  as  the 
ending  point  of  the  first  storm.  Figure  6  shows  the  events  for  a  complex  storm.  Even  if  the  temperature 
and  density  are  required  at  some  point  during  the  second  storm  it  is  important  to  start  the  temperature 
integration  at  the  commencement  of  the  first  storm  and  carry  it  through  into  the  second  storm,  since  the 
thermosphere  would  already  be  heightened  when  the  second  disturbance  began. 


2004  Storm  Geomagnetic  index  Dst 


Figure  6.  A  multiple  storm  during  2004,  showing  the  different  storm  events. 
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B.  Dst  Temperature  Equations 


Wilson  et  al.25  suggested  that  on  a  global  scale  the  storm-time  thermosphere  acts  like  a  large 
thermodynamic  system  that  never  strays  far  from  equilibrium.  From  an  analysis  of  GRACE  density 
measurements  Burke  et  al.24  further  argued  that  the  energy  input  to  the  thermosphere  can  be  treated  as  a 
large  driven-but-dissipative  thermodynamic  system,  which  can  be  described  by  differential  equations 
similar  to  that  of  a  resistor-inductor  circuit.  The  driver  is  the  magnetospheric  electric  field.  They  also 
demonstrated  that  Dst  and  storm-time  changes  of  the  exospheric  temperature  dTc  share  the  same  driver  but 
have  different  relaxation  time  constants.  By  eliminating  the  electric  field  term  from  the  two  equations 
Burke  et  al.24  established  the  following  relation  to  determine  exospheric  temperature  responses  as  a 
function  of  Dst. 


1 

f  o 

(1—  )dTco+S 

Dst,- 

1— 

Dsto 

T, 

l  TJ 

(8) 


where  ii  and  x2  represent  the  temperature  and  Dst  relaxation  e-fold  times.  From  an  analysis  of  the  GRACE 
data  during  2004  storms  Burke  et  al.24  obtained  values  of  Xi  =  6.5  hours  and  x2  =  7.7  hours;  the  slope  S  « 
1.58.  Using  these  values  in  Equation  (8)  results  in: 

dTc,  =  0.846  dTc0+ 1.58  [Dst,-  0.870  DstJ  (9) 

The  subscript  1  above  represents  the  value  at  one  time  step  beyond  the  subscript  0  point.  The  storm  event 
algorithm  described  above  was  used  to  determine  the  event  times  for  the  storm  being  considered.  Equation 
(9)  was  then  integrated  from  storm  commencement  time  until  storm  end  time,  producing  exospheric 
temperature  change  values  every  hour  throughout  the  storm  period.  These  temperature  change  values  were 
input  into  the  JB2008  model  to  represent  the  geomagnetic  storm  effects  at  all  points  throughout  the  storm. 

Comparisons  of  orbit  averaged  density  values  were  obtained  using  results  from  Equation  (9)  and  the 
CF1AMP  and  GRACE  accelerometer  densities.  Since  it  had  been  shown  that  the  Dst  index  was 
proportional  to  “global”  thermospheric  variations  it  was  decided  to  use  orbit  averaged  values  for  all  the 
comparisons.  Using  Equation  (9)  produced  good  correlations  of  the  JB2008  model  density  with  the 
accelerometer  data,  but  it  was  noticed  that  the  model  and  data  deviations  became  greater  as  the  maximum 
storm  magnitude  decreased  among  all  storms.  It  was  decided  to  re-determine  the  value  of  the  slope  S  while 
accepting  Burke’s  values  of  both  the  relaxation  parameters.  An  optimization  study  during  the  JB2008 
development  determined  that  these  r,  and  r,  values  were  the  best  to  use.  This  was  done  for  a  number  of 

storms  varying  from  minor  to  major.  The  slope  value  for  each  storm  was  optimized  by  minimizing  the 
differences  of  the  JB2008  model  orbit  averaged  density  ratios  using  Equation  (9)  with  the  orbit  averaged 
accelerometer  ratios  during  the  main  phase  region.  The  newly  determined  slope  for  each  storm  was  then 
plotted  as  a  function  of  the  storm  Dst  minimum  value,  and  also  plotted  as  a  function  of  the  ADst 
(minimum-maximum)  value.  The  Dst  minimum  values  produced  the  least  scatter  of  the  data.  Figure  7  is  a 
plot  of  the  main  phase  slope  with  respect  to  the  Dst  minimum  value. 

Equation  (10)  represents  the  new  quadratic  function  for  S  as  a  function  of  the  Dst  minimum  (DstMiN  )  value. 
If  DstMIN  <  -450  then  S  =  -1.40. 

S  =  -1.5050 xlO415  (DstMIN)2  -  1.0604 xl0~°2  DstMIN  -3.20  (10) 
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Storm  Main  Phase  Slope  vs  Dst  Minimum 


Figure  7.  Optimized  storm  main  phase  slope  S  values  as  a  function  of  storm  Dst  minimum. 


Using  this  new  slope  quadratic  equation  produced  very  consistent  results  for  storms  of  all  magnitudes. 
However,  a  few  additional  adjustments  had  to  be  made  to  produce  even  better  results.  It  was  discovered 
that  starting  the  dTc  integration  at  dTc  =  0  for  the  storm  commencement  time  sometimes  resulted  in  large 
negative  temperature  changes  at  storm  start.  This  was  due  to  the  fact  that  the  thermosphere  was  already  at 
a  slightly  heightened  temperature  state.  Therefore,  it  was  decided  to  start  the  dTc  integration  with  a  value 
equal  to  a  temperature  change  obtained  from  Jacchia’s  1970  geomagnetic  storm  equation  using  the  3 -hour 
a p  value  (with  a  6.7  hour  lag  time)  at  the  start  time.  Further  analysis  of  all  the  storms  showed  that  this 
produced  better  residts  than  using  an  initial  zero  value.  A  second  adjustment  during  this  main  phase 
analysis  occurred  during  sub-storms  when  the  Dst  variations  became  positive.  The  density  values  did  not 
drop  as  expected.  In  fact  the  accelerometer  and  HASDM  density  changes  during  these  time  periods 
continued  to  increase  even  though  the  Dst  value  was  increasing  during  these  short  main  phase  time  periods. 
Additional  equations  were  developed  for  these  time  periods: 

dTc,  =dTc0-SFACS(Dst,-Dst0)  (11) 

where  the  best  factor  SFac  was  found  to  be  0.3  for  all  storms.  The  S  value  was  obtained  from  Equation 
(10).  Since  S  is  negative  and  ADst  is  positive  during  these  time  periods  Equation  (11)  has  the  effect  of 
continuing  to  increase  the  temperature  change  (and  therefore  the  density)  even  though  Dst  is  increasing 
during  these  times.  Using  Equation  (11)  in  the  JB2008  model  produced  better  correlations  with  the 
accelerometer  data.  Finally,  it  was  noticed  that  a  small  lag  time  was  needed  to  better  represent  the  main 
phase  density  increase,  especially  during  small  storm  events.  It  was  determined  that  for  large  storms  (Dst  < 
-350),  moderate  storms  (-350  <  Dst  <  -250),  and  minor  storms  (-250  <  Dst)  lag  times  of  0,  1,  and  2  hours 
respectively  better  represented  the  main  phase  density  changes. 

The  recovery  phase  was  addressed  after  the  main  phase  equations  had  been  developed.  Equation  (9)  and 
(10)  were  initially  used  to  represent  the  recovery  phase  changes.  This  did  work  well  except  for  a  few 
outstanding  cases.  Each  storm  was  re-optimized  for  the  recovery  phase  by  optimizing  the  slope  for  this 
phase  only.  However,  the  recovery  phase  of  the  large  2003  multiple  storm  did  not  fit  the  accelerometer 
density  data  even  with  optimizing  Equation  (9)  just  for  the  recovery  phase.  It  was  decided  to  optimize 

r,  and  r  ,  for  this  phase.  After  many  trials  the  best  fit  for  the  2003  multiple  storm  was  r,  =  oo  and  r ,  =  I  . 
A  new  slope  was  then  obtained  for  this  storm,  and  the  resulting  equation  for  this  large  storm  was: 

14 

American  Institute  of  Aeronautics  and  Astronautics 


dTc,  =  1.00  dTc  +  0.13  Dst 


(12) 


The  next  step  was  to  determine  the  varying  slopes  for  storms  of  other  magnitudes.  Surprisingly  Equation 
(12)  was  found  to  be  the  best  representation  for  all  the  other  storms  representing  all  magnitudes.  This 
single  slope  value  was  excellent  for  the  entire  recovery  phase  up  to  the  recovery  slope  change. 

The  final  equation  fitting  was  for  the  period  covering  the  recovery  slope  change  to  the  end  of  the  storm.  It 
was  decided  to  use  the  simpler  Equation  (13)  below  since  it  was  supposed  that  for  this  time  period  the  ring 
current  had  disconnected  from  the  ionosphere,  which  meant  that  the  function  representing  the  ring  current 
energy  release  was  unknown. 

dTc,  =dTc0+S(Dst,-Dst0)  (13) 

The  resulting  slope  S  was  found  to  be  a  constant  -2.5  to  best  fit  all  the  storms.  It  was  found  sometimes  that 
dTc  became  negative  towards  the  “end”  time  of  the  storm  because  the  end  time  was  not  defined  correctly. 
To  make  sure  this  didn’t  occur  the  algorithm  sets  dTc  =  0  when  the  integration  step  produces  a  negative 
dTc.  Finally,  for  Dst  “non-storm”  periods  (Dst>-75)  JB2008  uses  Jacchia’s  1970  dTc  equation  as  a 
function  of  the  3-hour  ap  value.  When  the  JB2008  storm  computation  algorithm  has  determined  that  no  Dst 
storm  is  present,  then  if  ap  >  50,  a  value  of  50  is  used  for  the  dTc.  This  avoids  large  spurious  density 
increases  due  to  high  ap  values  when  no  storm  really  exists. 

C.  Dst  Modeling  Results 

Using  these  equations  for  each  of  the  3  different  storm  phases  results  in  very  good  comparisons  of  the 
JB2008  density  values  with  the  accelerometer  and  HASDM  values.  Figures  8  and  9  below  are  examples  of 
plots  of  model  density  ratios  during  two  major  storm  periods.  Yearly  average  density  values  were  obtained 
for  the  CF1AMP  and  GRACE  data.  The  displayed  CF1AMP  density  ratios  are  orbit  averaged  values  /  yearly 
average,  and  then  multiplied  by  1.17  to  adjust  to  the  F1ASDM  values.  The  17%  factor  is  based  on 
averaging  the  CF1AMP/F1ASDM  ratios  over  the  2001-2005  time  period.  A  factor  of  0.74  was  obtained  for 
the  GRACE/F1ASDM  ratios  based  on  all  data  from  2002  through  2005.  The  F1ASDM  values  plus  other 
model  values  are  orbit  averaged  (along  the  CF1AMP  or  GRACE  orbit),  with  all  ratios  based  on  each  year’s 
CF1AMP  (or  GRACE)  average  density  value.  Figure  8  shows  the  2004  major  storm  period  when  the 
GRACE  accelerometer  data  was  available,  and  Figure  9  shows  the  2003  major  storm  period  when  the 
CF1AMP  accelerometer  data  was  available.  The  F1ASDM  ratios  agree  extremely  well  with  the 
accelerometer  data  following  the  single  calibration  for  each  data  set.  The  JB2008  model  also  is  very 
consistent  with  the  density  changes  throughout  each  storm,  indicating  that  the  JB2008  model  temperature 
equations  are  working  extremely  well  for  these  orbit  altitudes  of  400  to  500  km.  The  MSIS  (NRLMSIS) 
density  values  are  mostly  low  at  storm  peak  times  during  the  largest  storms,  which  is  consistent  with  the 
results  previously  reported  by  Burke  et  al24.  The  Jacchia  70  (J70)  values  are  extremely  high  at  peak  storm 
times  because  they  are  based  on  single  ap  values  which  are  maxed  out  at  a  value  of  400  when  the 
magnetometers  are  saturated.  For  the  2003  storms  in  Figure  9  both  the  MSIS  and  J70  values  before  and 
after  the  storm  periods  are  much  too  high,  a  result  of  not  correctly  modeling  the  solar  EUV  during  this 
period  when  the  27-day  F 10  values  were  exceptionally  high. 

Finally,  Figure  10  shows  1 -standard  deviation  model  density  errors  as  a  function  of  storm  magnitude.  The 
values  were  obtained  as  percent  density  differences  from  the  calibrated  orbit  averaged  accelerometer  data, 
from  both  CF1AMP  and  GRACE,  and  the  different  model  orbit  averaged  values.  The  results  show  that  the 
JB2008  model  is  a  major  improvement  over  modeling  density  changes  during  large  geomagnetic  storms. 
The  F1ASDM  modeling  is  the  best  at  under  a  10%  sigma,  which  is  expected  since  it  accounts  for  real  time 
density  changes.  The  J70  modeling  is  the  worst  since  it  is  based  on  computing  a  density  from  a  single  3- 
hour  ap  value,  while  the  MSIS  model  uses  a  history  of  ap  values  for  57  hours  prior  to  the  time  of  interest. 
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2004  Dst  with  Orbit  Averaged  Density  Ratios 


Day  of  Year 


Figure  8.  Major  2004  storms  with  Dst,  ap  (left  scale) ,  and  density  (Rho)  ratios  displayed.  The  density 
ratios  are  based  on  orbit  averaged  model  density  values  /  GRACE  2004  density  average.  See  text  for 
additional  information. 


2003  Dst  with  Orbit  Averaged  Density  Ratios 


Figure  9.  Major  2003  storms  with  Dst,  ap  (left  scale),  and  density  (Rho)  ratios  displayed.  The  density 
ratios  are  based  on  orbit  averaged  model  density  values  /  CHAMP  2003  density  average.  See  text  for 
additional  information. 
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Orbit  Averaged  Model  Density  Errors 
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Figure  10.  Model  density  1-standard  deviation  errors  as  a  function  of  a p  ranges  representing  storm 
magnitudes.  Values  are  based  on  orbit  averaged  percent  density  differences  between  the  calibrated 
accelerometer  data,  from  both  CHAMP  and  GRACE,  and  the  different  model  values.  JB2006  uses  the 
same  geomagnetic  storm  modeling  as  J70. 


VII.  JB2008  Model  Usage 

A  detailed  description  of  the  model,  Fortran  source  code,  indices,  and  published  papers  can  be  obtained  at 
the  web  site  http://sol.spacenvironment.net/~ib2008/.  On  the  main  web  page  are  selections  for  an 
introduction  to  the  model,  all  the  papers  published  with  regard  to  the  JB2006  and  JB2008  development,  a 
separate  solar  and  geomagnetic  index  page,  the  Fortran  source  code  for  the  JB2008  model  and  driver 
program,  plus  selections  for  author  contacts,  figures,  and  the  SET  web  site.  On  the  indices  page  is  listed  a 
selection  for  the  Fortran  source  code  use  in  developing  the  new  temperature  index  from  the  Dst  values,  the 
new  temperature  index,  and  the  Dst  and  ap  indices.  A  short  description  of  the  Fortran  programs  follows. 

The  Fortran  source  code  for  the  JB2008  model  consists  of  two  source  code  files.  One  file  is  for  the  JB2008 
model  subroutine,  including  all  associated  subroutines  called  by  the  model.  The  inputs  and  outputs,  along 
with  the  publication  equation  references,  are  listed  within  the  code.  The  second  Fortran  source  file  is  a  test 
driver  program,  JB2008DRV,  with  associated  subroutines  required  to  call  the  JB2008  model  program. 
Also  included  are  a  solar  index  file,  the  geomagnetic  storm  temperature  index  file,  and  satellite  geographic 
position  test  input  and  output  files.  In  this  driver  program  the  solar  indices  are  obtained  from  the 
SOLFSMY  subroutine  that  uses  different  lag  times  for  retrieving  the  different  solar  indices.  A  test  input 
file  of  satellite  positions  is  used  to  compute  the  position  parameters  required  for  input  to  the  JB2008  model. 
Please  note  that  the  satellite  longitude  is  converted  to  a  right  ascension  value  required  by  JB2008.  A 
subroutine  for  computing  the  sun’s  position  is  also  included.  The  geomagnetic  storm  index  used  by 
JB2008  is  the  exospheric  temperature  change.  This  index  is  described  in  the  following  paragraph.  For 
every  input  test  satellite  position  the  JB2008  model  is  called  to  compute  the  density  value,  and  a  test  output 
file  is  written  with  the  density  and  indices  used  in  the  generation  of  the  density  value.  This  output  should 
be  compared  with  the  sample  output  provided  with  the  code  to  insure  the  program  is  running  correctly. 

The  geomagnetic  storm  temperature  index  used  by  JB2008  reflects  the  change  in  the  exospheric 
temperature.  This  temperature  change  is  computed  from  the  change  in  the  Dst  index  during  a  storm. 
Outside  of  storm  periods  this  index  is  computed  from  Jacchia’s  1970  equation  using  the  3-hour  ap  value. 
Since  the  temperature  change  is  global  in  nature  it  only  needs  to  be  computed  once  independent  of  any 
geography  (latitude,  longitude,  altitude)  when  new  Dst  values  are  obtained.  The  web  site  will  provide  the 
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continuously  updated  temperature  index  on  a  daily  basis  as  new  Dst  values  are  obtained  and  added  to  the 
Dst  fde.  The  temperature  index  file  is  then  used  as  an  input  to  the  JB2008  model.  Refer  to  the  test  driver 
JB2008DRV  source  code  for  its  usage.  Fortran  source  code  used  to  compute  this  temperature  change  is 
included  in  the  DTCMAKEDR  file.  The  file  has  a  driver  program  that  inputs  Dst  values,  determines 
geomagnetic  storm  event  times,  and  integrates  the  temperature  equations  as  a  function  of  the  Dst  changes 
throughout  an  entire  storm  period.  Refer  to  the  discussion  in  the  JB2008  model  publication  paper.  This 
computer  program  is  provided  for  reference  only  since  the  temperature  index  file  will  be  generated  on  a 
daily  basis  and  provided  at  the  web  site.  To  run  this  program  requires  the  Dst  file  and  the  3-hour  ap  file, 
both  of  which  are  also  found  on  the  web  site. 


VIII.  Conclusions 

The  following  results  have  been  obtained  using  the  new  JB2008  model: 

1.  Use  of  new  global  exospheric  temperature  equations  based  on  EUV  and  FUV  solar  indices  significantly 
improves  density  modeling,  especially  at  solar  minimum  times. 

2.  Density  standard  deviation  errors  during  non-storm  periods  have  been  reduced  by  over  5%  from 
previous  Jacchia  70  and  NRLMSIS  models. 

3.  Use  of  new  semiannual  density  variation  equations  using  multiple  81-day  averaged  solar  indices  now 
accounts  for  major  yearly  semiannual  density  changes  due  to  changing  long  term  EUV  heating. 

4.  Use  of  the  Dst  index  as  a  replacement  for  ap  greatly  reduces  density  errors,  especially  during  major 
geomagnetic  storm  periods.  This  error  reduction  is  from  over  60%  for  Jacchia  70  and  over  35%  for 
NRLMSIS,  to  16%  for  JB2008  during  major  storms. 

Significant  improvements  in  empirical  density  modeling  have  been  obtained  using  the  new  JB2008  model 
incorporating  new  solar  indices,  a  new  semiannual  variation  equation,  and  a  new  geomagnetic  index.  The 
new  model,  Jacchia-Bowman  2008  (JB2008)  provides  standard  deviations  of  approximately  9-10%  at  400 
km,  a  significant  decrease  from  16%  previously  obtained  using  the  Jacchia  70  model. 
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